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£f) , Abstract. We propose a first example of a simple classical field theory with 

nonholonomic constraints. Our model is a straightforward modification of a Cosserat 
rod. Based on a mechanical analogy, we argue that the constraint forces should 

^) | be modeled in a special way, and we show how such a procedure can be naturally 

implemented in the framework of geometric field theory. Finally, we derive the 
equations of motion and we propose a geometric integration scheme for the dynamics 
of a simplified model. 
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1. Introduction 

Motivation Nonholonomic constraints in mechanical systems have a long and illustrious 
history going back to the work of Hertz at the end of the nineteenth century. For 
nonholonomic constraints in field theories, the case is not so clear, and to the best of 
our knowledge, no convincing example has been proposed to this day. In this paper, 
we intend to rectify this omission by giving a simple example of a continuum theory 
with nonholonomic constraints. The basic model is that of a Cosserat rod, a special 
kind of continuum theory. This rod moves in a horizontal plane which is supposed to 
be sufficiently rough, so that the rod rolls without sliding. 

In the Cosserat theory one assumes that the laminae at right angles to the centerline 
are rigid discs. The nonholonomic Cosserat rod therefore touches the plane along a 
curve (rather than in an open subset of the plane, which would be the case for a fully 
three-dimensional continuum) and the nonholonomic constraint translates to the fact 
that the instantaneous velocity of these contact points is zero. 

Our example can also be modeled as the continuum limit of a nonholonomic mechanical 
system. Consider N rigid discs rolling vertically without sliding on a horizontal plane, 
and assume that these discs are interconnected by flexible beams of length £/N, as in 
figure [U Now let the number of discs go to infinity, while keeping the total length £ 
fixed: the result is the nonholonomic Cosserat rod. 

This mechanical model is interesting for a number of reasons. First of all, the 
nonholonomic field equations are derived by varying the action with respect to admissible 
variations, and this obviously requires the specification of a bundle of admissible 
variations, or equivalently, a bundle of reaction forces. In mechanics, this is commonly 
done by taking recourse to the principle of d'Alembert, which states that the virtual 
work of the reaction forces is zero. In field theory, this principle can be interpreted in 
a number of non-equivalent ways, and it is the mechanical model which will eventually 
determine our choice. 

Secondly, our model is a counterexample to the often-held belief that constraints in 
classical field theories are necessarily vakonomic. In sections 11.11 and 11.21 these two 
aspects are treated more in detail. 

Plan of the paper After giving a quick overview of jet bundle theory, we derive the 
Euler-Lagrange equations in the presence of nonholonomic constraints. Our treatment 
relies on the fact that the space of independent variables is a product of space and time, 
and that the fields are sections of a trivial fibre bundle. This is the case, for instance, 
for nonrelativistic elasticity. These assumptions allow us to split the jet bundle in a part 
involving spatial derivatives, and a part involving derivatives with respect to time. Using 
this natural splitting, we propose a bundle of reaction forces, based on the mechanical 
model (to be outlined in section [L2|) . As we shall see, these reaction forces are very 
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Figure 1. Geometry of the constrained rod 



similar to the ones used in mechanics. 

The remainder of the paper is then devoted to the study of a specific example of a 
nonholonomic field theory. First, we give an outline of the theory of Cosserat rods 
moving in the plane, and we pay particular attention to aspects of symmetry. In 
the second part, we then outline a suitable class of nonholonomic constraints, and we 
derive the equations of motion. The paper concludes with a brief foray into the field 
of geometric integration, where a simple explicit algorithm for the integration of the 
nonholonomic dynamics is proposed. 

1.1. Relation with other approaches 

In a number of papers [U [2], Bibbona, Fatibene, and Francaviglia contrasted the 
vakonomic and the nonholonomic treatments for classical field theories, and showed 
that for relativistic hydrodynamics only the former gives correct results. Another typical 
example of a vakonomic constraint is the incompressibility constraint in nonrelativistic 
fluid dynamics, treated by Marsden et al. [3]. Many more can be found in Antman's 
book [1] and in the papers by Garcia et al. [5] . 

In contrast, our field theory arises as the continuum limit of the vertically rolling 
disc, a textbook example of a nonholonomic mechanical system. These nonholonomic 
constraints survive in the continuum limit and hence provide a very strong motivation 
for the study of nonholonomic techniques in field theories. 

In previous papers (see 0, [3, El E] ) various theoretical frameworks were established for 
the study of nonholonomic field theories. Our model fits into these descriptions, but 
involves a number of additional ingredients which cannot be derived from these theoretic 
considerations alone. In particular the bundle of reaction forces takes a special form, 
motivated by similar definitions from mechanics. 

It should also be noted that similar theories as ours were explored before by Vignolo and 
Bruno (see [10J). They considered constraints depending only on the time derivatives of 
the fields, and their resulting analysis is therefore more direct. However, the underlying 
philosophy is the same: the constraints are "(...) purely kinetic restrictions imposed 
separately on each point of the continuum" . 
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1.2. Modeling the constraint forces 

Nonholonomic mechanical systems The mechanical background is not essential for the 
description of the continuum theory, but rather serves as a justification for some of our 
definitions. In particular, it provides a number of valuable clues regarding the type 
of constraint forces needed to maintain such a nonholonomic constraint. Let S be the 
configuration space of the vertically rolling disc, so that the configuration space for 
the entire model, consisting of N discs, is the product space S N . Denote by tp% the 
constraints of rolling without sliding imposed on the ith wheel; if?^ is a function on 
TS N . 

With these conventions, a motion of the system is a curve t i— > c(t) in S , and a variation 
of such a motion c is then a vector field (X 1; X 2 , ■ ■ ■ , X N ) on S N along c, i.e. a collection 
of maps Xi : K — > TS such that Xi(t) G T&rqS for all i — 1, . . . , N, where c l := pr^ o c. 

Let us now consider the one- forms $?k := J*(dip?^), where J is the vertical 
endomorphism on TS N . In geometric mechanics, linear combinations of these one-forms 
represent the possible reaction forces at the ith wheel; the bundle F, defined as 

F:=<$f 1) >©<$f 2) >©---©<$f 7V) > 

then represents the totality of all reaction forces along the rod. In coordinates, the 
one-forms $?-, are given by 

dip a ( d(p a \ 

*«=i% d *' ( = ^si dyh ) ta « n '= i -- w - ^ 

Here, (yu),y(i)) is a coordinate system on the ith. factor of TS N . Note that there is no 
summation over the index i in flTJ), and that the summation over individual coordinates 
is implicit, as shown in the term between brackets. 

Knowing the precise form of the bundle of reaction forces F is important because the 
nonholonomic equations of motion are derived by varying the action with respect to 
admissible variations. Moreover, the principle of d'Alembert shows us that a variation 
is admissible if it belongs to the annihilator of F, i.e. a variation (X l5 . . . ,X^) of c is 
admissible if 

(Xi(t),a{c{t))) = Q, for all(i,i) G {l,...,iV} xR and a G F, (2) 

where Xj is a lift of Xi to T(TS) such that Tt$ ° Xi = Xi. Note that there is again no 
summation over i. In coordinates, this is equivalent to 

Vi^- for alii = 1,...,JV, (3) 



where we have written X,- = v,-. 



a 



9y(i) 
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In the next paragraph, we will let the number N go to infinity, while keeping the length P. 
constant. The result is a field theory, and a reaction force will be a continuous assignment 
of a one-form on TS to each point of the centerline of the rod. This definition will be 
the starting point for our treatment in the main body of the text; once we know the 
bundle of reaction forces, we can then derive the nonholonomic field equations. 

The continuum model In the continuum limit, a field is a map <j) from [0, £] x R to S. It 
is customary in classical field theory to view these fields as sections of a trivial bundle it, 
whose base space is [0, £] x R, and with standard fibre S. The role of the velocity space 
is then played by the first jet bundle J 1 tt, and the constraints (pry from the previous 
paragraph are replaced by a constraint function ip a on J 1 it. 

In field theory, a variation of a field is now a map X : [0,l]xl-> TS with the 
property that X(s,t) G T^^S; in other words, a vector field along 0. Taking our cue 
from (IHl), we say that a variation X is admissible if the following holds (in coordinates): 

dip a 



where we have written X(s,t) = X(s,t) a -^. 

This condition can be rewritten in intrinsic form by using the following observation: as 
we shall show below, there exists a natural isomorphism between the first jet bundle and 
the product bundle R x [^(M^) x s TS}. Now, let J be the vertical endomorphism 
on TS. This map has a trivial extension to the whole oflx [J l (M, S) ^sTS], and by 
using the natural isomorphism with J l 7t, we obtain a map J* from T*(J 1 7r) to itself. 
The bundle F of constraint forces is then generated by the forms $ a := J*(df a ). The 
similarity with the mechanical case is obvious. 



2. Lagrangian field theories 

In this paper, we will mostly be concerned with the description of elastic bodies. The 
geometric description of these theories is well known and we refer to [Til [T2] for more 
information. 

Let M be a smooth n-dimensional compact manifold, and let 5* be a general smooth 
m-dimensional manifold, with n < m. The points of M are "material points", labelling 
the points of the body, whereas S is the physical space in which the body moves. In 
most cases, S will be the Euclidian space M. 3 , whereas M can be one-, two-, or three- 
dimensional, corresponding to models of rods, shells, and three-dimensional continua. 
Furthermore, M is assumed to be oriented, with volume form t]m- 

On M we consider a coordinate system (x l ) , i = 1, . . . , n, such that r\ M can locally be 
written as tjm = dx 1 A- ■ ■ Adx n , and on S we take a coordinate system (y a ), a — 1, ... ,m. 
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2.1. The bundle picture 

In this section, we give a brief overview of the theory of jet bundles. For an introduction 
to classical field theory using jet bundles, we refer to [T3l [HI US] and the references 
therein. 

Consider the fibre bundle it : Y — > X, where X :=IRxM, Y := X x S, and the 
projection ir is the projection onto the first factor. The manifold X is equipped with a 
coordinate system denoted by (a; M ), where fi = 0, . . . , n, and such that x° := t. Similarly, 
S has a coordinate system (x^, y a ) which is adapted to the projection in the sense that 
the projection tc is locally given by pr 2 : (x M , y a ) h-» (x m ). Note that M is equipped with 
a volume form 77 = dt A 77a/, which we write in coordinates as d n+1 x := dx° A • • • A dx n . 
We will employ the following short-hand notation: 

d% := -^J 77 = (-lfdxo A • • • A dx^ 1 A dx^ +1 A • • • A dx n . 

The first jet bundle </% is the appropriate stage for Lagrangian first-order field theories. 
Its elements are equivalence classes of sections of 7T, where two sections are said to be 
equivalent at a point x G X if they have the same value at x and if their first-order 
Taylor expansions at x agree. The equivalence class of a section at x is denoted by 
jl4>. Hence, J x n is naturally equipped with a projection 7ri j0 : J 1 it — ► Y, defined by 
Trifl{jl4>) = </>(x), and a projection n 1 : J 1 ^ — > X, defined by 7Ti(j^0) = x. The induced 
coordinate system on </% is written as (x' I ,y a ;yf l ). 

Furthermore, we define the manifold J l (M,S) of jets of mappings from M to S as 
the first jet manifold of the trivial bundle pr x : M x S — > M. The usual jet bundle 
projections Hi and 7r 10 induce projections n M , onto M, and n s , onto S, respectively. 

Because of the special structure of the bundle n, viz. the fact that X is the product 
RxM and that it is trivial, J lr K can be written in a special form. Recall that the fibre 
coordinates of J l ir represent the derivatives of the fields with respect to space and time: 
the decomposition of lemma [27T1 then provides an invariant way of making the distinction 
between time derivatives and spatial derivatives. Such an invariant decomposition is not 
possible for general jet bundles. 

The bundle IR x [J^-^M^) x s TS], a fibered product, consists of triples (t, K, v) such 
that its (k) = t(v), where r : TS — > S is the tangent bundle projection. It is 
equipped with a projection tx onto Y defined as n(t, k, v) = (t, 7t m (k); t{v )). Moreover, 
n : R x [J\M, S)x s TS]^Y is an affine bundle. 

Lemma 2.1. The first jet bundle J l n is isomorphic, as an affine bundle over Y = 
RxMxS,toRx [J l {M, S) x s TS]. 

Proof: A more general statement can be found in section 6B of |16j . 

Take any point (t, m,s) mWx M x S and consider a 1-jet 7 such that vr 1)0 (7) = (t, m, s). 
An alternative interpretation of 7 is that of a linear map 7 : T(^ m )(R x M) — > T S S. 
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Consider now the map \P(t )ms ), mapping 7 to the element of R x [^(M, S) x s TS] given 
by 



*./.„,,,(;-) ■ (/.-(O,.-).- ( — 



; "m 



where m and t are the zero vectors in T m M and in TJJR, respectively. It is easy to 
check that the map \I/ is an isomorphism of affine bundles. □ 

In coordinates (y a , y a ) on TS and (x l , y a ; y") on J 1 (M, S), the isomorphism of lemma [2TT1 
is given by (a^, y a ; y°) \-> (*; z*, ?/ a , j/JJ = y 4 a ; ?/ a , y a = j/g). 

For future reference, we remark here that the tangent bundle TS is equipped with a 1-1 
tensor field, denoted by J : T(TS) —> T(TS), and given in coordinates by 

An intrinsic definition can be found in [T7]. The adjoint of J will be denoted by J* and 
is a map from T*(TS) to itself defined by (J* (a), v) = (a, J(v)), for all v G T(TS). 

In section |3j we will encounter higher-order field theories, in particular of order 2. In 
order to be able to deal with this type of field theory, we introduce the manifold J k 7f of 
kth order jets. The elements of J h n are again equivalence classes of sections of 7r, where 
two sections are equivalent at x G X if they have the same value at x and if their Taylor 
expansions at x agree up to the kth order. The kth order jet bundle is equipped with 
a number of projections itk,i '■ J k ^ — * J ft (where I < k), constructed by "truncating" 
to order I the Taylor expansion defining an element of J k ir. A detailed account of jet 
bundles is provided in [IS] . 

As a matter of fact, we will only need the third-order jet manifold J 3 n. A natural 
coordinate system on J 3 7r is given by (x^,y a ; y"; y^ v ; y^ UK ), for a = l,...,m and 
[i,u,K = Q,...,n, with the convention that 

Vfiv = Vv^ an d V^uk = Va{^vn) 

for any permutation o of the three indices (expressing the commutativity of partial 
derivatives). 



2.2. Covariant field theories of first and second order 

2.2.1. First-order field theories The geometry of first-order field theories has been 
studied by many authors (see [3j HH [151 EH HH] and the references therein for a non- 
exhaustive survey) and is by now well established. In this section, we recall some basic 
constructions. 

Consider a first-order Lagrangian L : J lr K — > M.. There exists an (n + l)-form Gl on 
J 1 TV, called the Cartan form. Different intrinsic constructions of 0£ can be found in 
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[131 El IB] , but here its coordinate expression will suffice: 

L = Ld n+1 x + — (dy a - y a v dx v ) A d%. (4) 

Let us also define the (n + 2)-form Ql '■= — d0L on J lr K. 
Let S be the action functional defined as 

S(0) = f LO'Vto, (5) 

for each section of 7r with compact support U. We now look for critical points of this 
functional under arbitrary variations, which are defined as follows. 

Definition 2.2. An infinitesimal variation of a field <fi defined on U is a vertical vector 
field V defined in a neighbourhood of<p(U) in Y , with the added restriction that V(y) = 
for all y E (f>(dU). 

An infinitesimal variation gives rise to a local one-parameter group of diffeomorphisms 
$ e : Y — > Y defined in a neighbourhood of 4>{U). The fact that V vanishes on <p(dU) 
implies that $ e (y) = y for all y G <p(dU). The composition of $ e with a section of 7r 
is hence a new section of it, denoted by </) e . 

The critical points of S therefore satisfy 



d 



dL d dL 



= S 5 W >(»..«) _ o = J a [ W - ^ W J ™~,. (6) 

As the variations V are arbitrary, we obtain the usual Euler-Lagrange equations for a 
section 6 of tt: 



dL d dL 



dy a dx^ dy a 



(J 2 0) = 0. 



These partial differential equations can be rewritten in intrinsic form by means of the 
Cartan form (see [131 H [19]): 

Theorem 2.3. A section of i\ is a critical 'point of the action S, or, equivalently, 
satisfies the Euler-Lagrange equations, if and only if 

{]^y{i w n L ) = o (7) 

for all vector fields W on J Xr n. 

Symmetries and Noether's theorem Let G be a Lie group acting on X by 
diffeomorphisms, and on Y by bundle automorphisms. For g e G, consider the bundle 
automorphism $ g : Y — > Y with base map f g :X^X. The prolongation to J l ii of 
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the action of G is then defined in terms of bundle automorphisms j 1( £ 5 : J 1 ^ —* J 1 ^, 
defined by 

J 1 $ 9 (j^)=4 (x) ($ 9 o0o/- 1 ). 

Consider an element £ of g and denote the infinitesimal generator of the prolonged 
action corresponding to £ by O 1 ^- Note that (,ji n is just j 1 ^, the prolongation of 
the infinitesimal generator on Y corresponding to £. We recall that if £y is given in 
coordinates by 



dxf* ' dy a 
then its prolongation j £y is defined as 



<9a^ ' <9y a \da;^ u dx^ J dy 



dx^ ) dy® 

We say that a Lagrangian L is invariant under the prolonged action of G if Lo j 1 ^ = 
(det[-D/ g ]) _1 L for all p G G and 7 G J 1 ^, where det[D/ 9 ] is the Jacobian of f g . This 
condition can be expressed concisely as (j 1 ^ g )*(Lrj) = Lrj. It can be shown that 
invariance of the Lagrangian implies invariance of the Cartan (n + l)-form, expressed 
as j^^Ol = @l for all g G G, or, infinitesimally, 

ST^Bl = 0. (8) 

Let L be a G-invariant Lagrangian. Associated to this symmetry is a map defined as 
Jf := 1 ?rJ@£- We now define the momentum map J L : J 1 ^ — > ^"(J 1 ^) Cg) g* by 
(j L ,^ = Jt . The importance of the momentum map lies in the following theorem, 
which we have taken here from [T4"l thm. 4.7]: 

Proposition 2.4 (Noether). Let L be an invariant Lagrangian. For all £ G Q, the 
following conservation law holds: 

d[(j V)*41 = °> 
/or a// sections <f> of 7r that are solutions of the Euler- Lagrange equations |7|). 

A comprehensive account of symmetries in classical field theory can be found in [20] . 

2.2.2. Second- order field theories Many of the Lagrangians arising in elasticity are of 
higher order. In particular, we will encounter a second-order model in section [31 In some 
papers (see [HJ [21] and the references therein) a geometric framework for second-order 
field theories has been developed and we now recall a number of relevant results. 



A class of nonholonomic kinematic constraints in elasticity 



10 



A second-order Lagrangian is a function L on J 2 tt. The corresponding second-order 
Cartan (n + l)-form is a form on J 3 tt, whose coordinate expression reads 



e, 



+ 



dL 



dL 



L- 



dyl dx^ \dyl 



dL 



IVjX 



dy a A d n x v + j—dy a u A d% 
dy a v 



dy a u 



yl + 



dL 



dx^ \dy% 



Vl 



dy a y^ 



d n+l 



X. 



(9) 



Many results from the previous section on first-order field theories carry over 
immediately to the higher-order case. The action 5* is defined as 

S(<j>) = [ Ltfflr,, 

where is again a section of 7r with compact support U. A section is a critical point 
of this functional if and only if it satisfies the second-order Euler-Lagrange equations: 



dL 



dL 



+ 



dy a dx>* ydy^J dxMx^ \dy® 



dL 



y fiv / . 



O'V) = o. 



(10) 



There also exists an intrinsic formulation of the Euler-Lagrange equations. We quote 
from [21J: 

Proposition 2.5. Let L be a second- order Lagrangian. A section <ft of ir is a solution 
of the second-order Euler-Lagrange equations if and only if (j 3 4>)*(W1Ql) = for all 
vector fields W on J 3 ir. Here, Ql := — dO^ is the second-order Poincare- Cartan form. 

Remark 2.6. It should be noted that there always exists a Cartan form for higher- 
order field theories, but that uniqueness is not guaranteed (contrary to the first-order 
case). However, by imposing additional conditions, Saunders [18] was able to prove 
uniqueness for second-order field theories. This unique form, given in (Q, was derived 
by Kouranbaeva and Shkoller [2T] by means of a variational argument. o 



The action of a Lie group G acting on Y by bundle automorphisms gives rise to a 
prolonged action on J 2 n. If a Lagrangian is G-invariant with respect to this action, then 
the momentum map J L E fi n (J 3 7r) <g> g*, defined as (J L ,£) = J^, where J^ = 3 7r©L, 
gives rise to a conservation law: d[(j 3 0)* Jp] = for all sections (f> of 7i that are solutions 
of the Euler-Lagrange equations ffTUj) . 



2.3. Nonholonomic field theories 

2.3.1. The field equations We now derive the Euler-Lagrange equations in the presence 
of nonholonomic constraints. The nonholonomic problem involves the specification of 
two distinct elements: the constraint manifold C, and the bundle of reaction forces F. 
Following Marie [22J , we impose no a priori relation between C and F. 

The constraint manifold C is a submanifold of J l n of codimension k and represents 
the external constraints imposed on the system. For the sake of definiteness, we will 
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assume that C projects onto the whole of Y (i.e. vr 10 (C) = Y), and that the restriction 
(7Tio)ic '■ C —>■ Y is a fibre bundle. This need not be an affine subbundle of J l 7t. For 
the benefit of clarity, C is assumed to be given here by the vanishing of k functionally 
independent functions <f a on J lr K\ 

C := {7 G JV : /( 7 ) = for a = 1,..., A;}. 

The treatment can be easily extended to the case where the (p a s are only locally defined. 

Secondly, we assume the existence of a fc-dimensional codistribution F on J 1 ^, along 
C, of reaction forces. The elements of F are maps a : C — > T*S such that a (7) G T*S, 
where s = (pr 2 o ^1,0) (7)- If we denote by tits '■ J 1 ^ — ► TS the composition 
kts '■= P r 3 ° ^? where \1/ is the isomorphism defined in lemma l2?lj then the elements 
of F can equivalently be viewed as one-forms along the projection ttts- By pull-back, 
these one-forms then induce proper one-forms defined along C. In local coordinates, an 
element a of F can be represented as a = A a (x^,y a ,y^)dy a , where the A a are local 
functions on C. 

We define the annihilator F° of F as the following subbundle of TJ 1 ^ along C: for all 

7 eC, 

F° := {f 7 G TyJ 1 ^ : (a 7 , t> 7 ) = for all a 7 e F 7 }. 
An arbitrary element v 7 of F° has the following form: 

Here, we have chosen a basis of sections A"dy a of F . No further restrictions are imposed 
on the coefficients v^ and u", but note that this is not the end of the story; see the 
appendix. 

The local work done by a "force" a along a variation V (see definition 12.21) is then given 
by the pairing (0(7), j 1 V( , y)), where 7 G Imj 1 ^, and the global work done at time t by 
the integral JmU 1 ^)* (&,j 1 V) r lM, where 4>t is the instantaneous configuration defined 
by (j) t {u) := <p(t,u). 

Definition 2.7. A variation V of a field <f> defined over an open subset U with compact 
closure is admissible if (J 1 <f>)*(j 1 Vlcx) = for all a G F. 

Definition 2.8. A local section <p ofir, defined on an open subset U C X with compact 
closure, is a solution of the nonholonomic problem determined by L, C, and F if 
j l (p(U) C C and (TJJ) holds for all admissible variations V of (f). 

It follows from ([6]) that a local section <p is a solution of the nonholonomic problem if it 
satisfies the nonholonomic Euler-Lagrange equations: 



dL d dL 



dy a dx^ dy a 



(fcf>) = Aoi^O'V) and ^(j'V) = 0. (11) 
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Here, X a are unknown Lagrange multipliers, to be determined from the constraints. 
This is proved below. 

Theorem 2.9. Let <p be a section of it. Iflmj l <f) C C, then the following assertions are 
equivalent: 

(a) (p is a stationary point of the action (EJj under admissible variations; 

(b) cf) satisfies the Euler- Lagrange equations (TO)); 

(c) for all n-vertical vector fields V such that j}V G F° for all 7 G C, 

(jV)*0'Vj^) = 0. (12) 

Proof: Let us first prove the equivalence of (jaj) and (jcj). For arbitrary, not necessarily 
admissible variations, the following result holds (this is equation 3C.5 in |14J): 

%-Sfa) =- /(jV)WJ^). 
de ,=o J x 

For admissible variations, we have therefore 

X 

Now, we may multiply V by an arbitrary function on X and this result will still hold 
true. The fundamental lemma of the calculus of variations therefore shows that 

(jV)*(ji^J^) = 0, (13) 

for all admissible variations V. By using a partition of unity as in [14], it can then be 
shown that (!T2l holds for all 7r-vertical vector fields V such that (j 1 V,a) = for all 
a G F. 

The equivalence of (jb]) and (Jcj) is just a matter of writing out the definitions. In 
coordinates, the left-hand side of (TT21) reads 

;iX\*fAiifin.\ - I/O/-! a\ f dL (a\ x\_ d dL fAx\\A»+i, 



and this holds for all variations j 1 V G F°. Therefore, if (p satisfies (|12[) . then there exist 
k functions X a such that 

The converse is similar. □ 

Remark 2.10. In f fT2|) . only prolongations of vertical vector fields were considered, 
whereas in similar expressions in theorem [231 and proposition [231 arbitrary vector fields 
occurred. This is due to the fact, also mentioned in the appendix, that only vertical 
variations are considered. 
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For the derivation of the nonholonomic Euler-Lagrange equations, it is enough to 
consider only vertical variations. However, if one wants to prove the nonholonomic 
Noether theorem for symmetries that act nontrivially on the base space, one needs an 
expression like (Tl2|) but with j l V replaced by a vector field W which is not necessarily 
the prolongation of a vertical vector field. As we point out in the appendix, this can be 
done, but then one needs to modify the definition of F . o 

2.3.2. The Chetaev principle In section 12.3. H we defined reaction forces as certain 
maps from C to T*S, but their exact nature was left unspecified. We now conclude our 
derivation of the nonholonomic field equations by proposing a concrete definition for 
these reaction forces. As will become clear in a moment, this definition is formally 
identical to the one used in mechanics; roughly speaking, the reaction forces are 
constructed by composing dip a with the vertical endomorphism J on TS, which is what 
one might call the Chetaev principle. 

Indeed, the vertical endomorphism J on TS trivially extends to a (1, l)-tensor J on 
R x [J x (Af, S) x s TS\, defined as 

J(a,/3, 7 ):=(0,0,J( 7 )), (14) 

where a e T t *R, (3 G T*J\M,S), and 7 e T*{TS) (and where (t,u,v) e R x 
[J^M, S) x s TS]). We denote the adjoint of this map as J*. 



Let ip a be the k constraint functions. By means of the isomorphism \1/ of lemma I2.1[ 
these functions induce k functions on R x [.^{M, S) x s TS], which we also denote by 

Lf a . 

Definition 2.11. The bundle of reaction forces F is the co- distribution on J lr K locally 
generated by the following forms: F = Span($ Q ), where 

$ a : = ry*(d<^ Q )]. 

In local coordinates on J 1 ^, the generating forms $ a are given by 

^ = ^-dy a . (15) 



This corresponds to the coordinate expressions based on the mechanical analogue: 
compare, for instance, with (pQ). Again, we emphasize that there is an obvious distinction 
between spatial derivatives and derivatives with respect to time. 

Using these reaction forces, the nonholonomic Euler-Lagrange equations become 



dL d dL 



tfV) = AainrO'V) and ^ a (iV) = o, 



dy a dx 11 dy^ 

where the X a are again a set of unknown Lagrange multipliers, to be determined from 
the constraints. 
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Figure 2. Geometry of the constrained rod 



3. A Cosserat-type model 

The theory of Cosserat rods constitutes an approximation to the full three-dimensional 
theory of elastic deformations of rod-like bodies. Originally conceived at the beginning 
of the twentieth century by the Cosserat brothers, it laid dormant for more than fifty 
years until it was revived by the pioneers of rational mechanics (see [HJ §98] for an 
overview of its history). It is now an important part of modern nonlinear elasticity and 
its developments are treated in great detail for instance in [I], which we follow here. 

A Cosserat rod can be visualised as specified by a curve s i— > r(s) in M 3 , called the 
centerline, to which is attached a frame {di(s), d 2 (s), d 3 (s)}, called the director frame 
(models with different numbers of directors are also possible). The rough idea is that 
the centerline characterizes the configuration of the rod when its thickness is neglected, 
whereas the directors model the configuration of the laminae transverse to the centerline. 
In the Cosserat theory, the laminae are assumed to deform homogeneously, and therefore 
the specification of a director frame in M 3 , fixed to a lamina, completely specifies the 
configuration of that lamina. 

In the special case where the laminae are rigid discs at right angles to the centerline, 
one can choose the director frame {d l5 d 2 ,d 3 } to be orthogonal with, in addition, d 2 
and d 3 of unit length (attached to the laminae) and di aligned with r'(s), the tangent 
vector to the centerline. If, in addition, the centerline is assumed to be inextensible, so 
that we may choose the parameter s to be arclength, di is also of unit length and the 
director frame is orthonormal. In this case, the specification of, say, d 2 is enough to 
determine a director frame: putting d x = r', we then know that d 3 = d x x d 2 . Here 
and in the following, a prime (') denotes derivation with respect to s. 

Here, we will consider the case of a Cosserat rod with an inextensible centerline and 
rigid laminae. In addition, we will assume that the centerline is planar in the (x, y)- 
plane, which will allow us to eliminate the director frame almost completely. The result 
is a Lagrangian field theory of second order, to which the results of section 12.21 can be 
applied. 
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3.1. The planar Cosserat rod 

Consider an inextensible Cosserat rod of length £ equipped with three directors. If we 
denote the centerline at time i ass ^ r(t,s), inextensibility allows us to assume that 
the parameter s is the arclength. Secondly, we can take the director frame {di,d2,d3} 
to be orthonormal, such that di is the unit tangent vector r'. We will not take the effect 
of gravity into account. 

In addition, we now assume that the centerline is a planar curve moving in the horizontal 
(x, y)-plane, i.e. r(t,s) can be written as (x(t,s),y(t,s),0). We introduce the slope 
<p(t,s) of the centerline as (cos </?, sin ip) = (x'(t,s),y'(t,s)). Furthermore, we define the 
angle 6(t,s), referred to as the torsion of the rod, as the angle subtended between e z 
and d 3 . The director frame is completely determined once we know the slope tp(s,t) 
and the torsion 9(s,t). 

The specific constraints imposed on our rod model therefore allow us to eliminate the 
director frame in favour of the slope p> and the torsion 6. Furthermore, as we shall see, 
the slope tp is related to the curvature of the centerline. Note that, in formulating the 
dynamics, we still have to impose the inextensibility condition (x') 2 + (y') 2 = 1. 

Remark 3.1. Note that 6 has nothing to do with the usual geometric concept of torsion 
of a curve in M 3 , and neither is 9 related to the concept of shear in (for example) the 
theory of the Timoshenko beam. o 

3.2. The dynamics 

As the director frame is orthonormal, there exists a vector u, defined by (!■ = u x dj, 
called the strain or Darboux vector. With the conventions from the previous section, u 
takes the following form: 

u = 0'di + <p'e z . 

(u can be thought of as an "angular momentum" vector, but with time-derivatives 
replaced by derivatives with respect to arclength.) 

The dynamics of our rod model can be derived from a variational principle. The kinetic 
energy is given by 



T = \J (p(s)(x 2 + y 2 ) + a6 2 )ds, 



where a is an appropriately chosen constant. Here, the mass density is denoted by p, 
and will be assumed constant from now on. 

For a hyperelastic rod, the potential energy is of the form 

r-C 



V = W(u 1 ,u 2 ,u 3 )ds, 
Jo 
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where W{u\, 1X2,11%) is called the stored energy density, and the Ui are the components 
of u relative to the director frame: Ui — u • dj. In the simplest case, of linear elasticity, 
W is a quadratic function of the strains: 

W( Ul , u 2 , u 3 ) = X - [K x u\ + K 2 u\ + K 3 u 2 3 ) . (16) 

We will not dwell on the physical interpretation of the constants Ki any further (in 
this case, they are related to the moments of inertia of the laminae). If the rod is 
transversely isotropic, i.e. if the laminae are invariant under rotations around di, we 
may take K 2 = K3. The potential energy then becomes 

V = lf(P(0') 2 + KK 2 )ds, (17) 

where k is the curvature of the centerline, i.e. k 2 = (f') 2 = (x") 2 + {y") 2 , and where 
we have put (3 := K\ and K := K 2 . Models with a similar potential energy abound 
throughout the literature and are generally referred to as the Euler elastica. For more 
information, see 1 23 1 and the references therein. 



3.3. The second-order model 

Having eliminated the derivative of the slope ip from the stored energy density, we end up 
with a model in which the fields are the coordinates of the centerline (x(t, s),y(t, s)) and 
the torsion angle 8(t, s). This model fits into the framework developed in section I2.2.2J 
the base space X is M. 2 , with coordinates (t, s) and the total space Y is X x IR 2 x S 1 , 
with fibre coordinates (x,y,8). 

The total Lagrangian now consists of the density of kinetic energy minus that of potential 
energy, as well as an additional term enforcing the constraint of inextensibility, and can 
be written as 

L = \ (x 2 + ij 2 ) + ^0 2 -\ Wf + Kn 2 ) - \ V {{x'f + {y') 2 - 1) , (18) 

where p is a Lagrange multiplier associated to the constraint of inextensibility. The field 
equations associated to this Lagrangian take the following form: 

px + Kx"" = §- 8 {j)x') 

py + Ky"" = l(py>) (19) 

a'e-pe" = 0, 

to be supplemented with the inextensibility constraint 

{x'f + {y') 2 = 1, (20) 

which allows to determine the multiplier p. Note in passing that the dynamics of the 
centerline and the torsion angle 9 are completely uncoupled. This will change once we 
add nonholonomic constraints. 
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3.4- Field equations and symmetries 

We recall the expression (Q for the second-order Cartan form. If a Lie group G is acting 
on Y by bundle automorphisms, and on J 3 ir by prolonged bundle automorphisms, there 
is a Lagrangian momentum map jf = G 3 7r-I^L> as described in section I2T21 We now 
turn to a brief overview of the symmetries associated to the rod model introduced in 
the previous section. For an overview of symmetries in the general theory of Cosserat 
rods, see [2] 



Translations in time The Lie group 1R acts on X by translations in time: $ e : {s, t) h- > 
(s, t + e). The Lagrangian is invariant and the pullback to X (by a solution j 3 of the 
field equations) of the momentum map associated to the infinitesimal generator J^ is 
given by 

(f(f))*j[ = \(px' - Kx'")x + (py' - Ky'")y + f39'9 + K(x"x' + y"y')] dt (21) 

+ [ P -{x 2 + y 2 ) + ^9 2 + f «z") 2 + {y") 2 ) + §(*0 a + l^'f + (V'f - 1)]^ 

^^ v ' 

e 

where we have introduced the energy density S. By taking the exterior derivative of 
( |2T|) and integrating the conservation law d[(j 3 0)*Jf] = over [0, £] x [t ,ti] C M. 2 , we 
obtain 

E{tx) - E(t ) = I J \{px' - Kx'")x + (py' - Xy w )y + P°'0 + K(x"x + y"y)] * dt, (22) 



/o 



J o 



where E(t) = f Q £ds is the total energy, which is conserved if suitable boundary 
conditions are imposed. This is the case, for instance, for periodic boundary conditions 
or when both ends of the rod can move freely, i.e. when 

px' - Kx'" = py' - Ky'" = and x" = y" = 9' = at s = 0, L 

Spatial translations Consider the Abelian group M? acting on the total space Y by 
translation, i.e. for each (a, b) € M 2 we consider the map $( a ,&) : (s, t; x, y, 6) \— > 
(s,t;x + a,y + b,9). The Lagrangian density is invariant under this action and the 
associated momentum map is 

{f4>)*J(vi,v2) = -p(vii + v 2 y)ds - (v t px' - v t Kx'" + v 2 py' - v 2 Ky'")dt 

for all (fi,f2) G M. 2 . Again, under suitable boundary conditions, Jh , gives rise to a 
conserved quantity, namely the total linear momentum of the rod. 

Similarly, S 1 acts on Y by translations in 9, with infinitesimal generator of the form J|, 
and the corresponding momentum map is 

(j 3 0)* jf = -(39'dt - a9ds. 
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The ensuing conservation law is given by ot9 = j39" and, hence, is just the equation of 
motion for 9. 



Spatial rotations Finally, we note that the rotation group SO (2) acts on Y by rotations 
in the (x, j/)-plane. The infinitesimal generator corresponding to 1 G so (2) = R is given 



by y-^; — x-§-\ its prolongation to J 3 n is 

d d d d . d , d 

G 3 7T = Vtt- ~ x— + y— - x— + yT n -x 7 — + ---, 
ox ay ox ay ox' ay 

where the dots represent terms involving higher-order derivatives. As 0^ is semibasic 
with respect to 7r 31 , these terms make no contribution to the momentum map. The 
momentum map is given by 

(j 3 0)* Jf = [-x(-py' + Ky") + y(-px' + Kx") - K(x"y' + y"x')} dt + p(xy - yx)ds, 

leading to the conservation of total angular momentum. Note that the angular 
momentum does not involve 8, in contrast to the corresponding expression in more 
general treatments of Cosserat media. This is a consequence of the fact that we defined 
the action of 50(2) on Y to act trivially on the 9 part. 



3.5. A nonholonomic model 

Consider again a Cosserat rod as in the previous section. The constraint that we are 
now about to introduce is a generalization of the familiar concept of rolling without 
sliding in mechanics: we assume that the rod is placed on a horizontal plane, which we 
take to be perfectly rough, so that each of the laminae rolls without sliding. 

However, as the Cosserat rod is also supposed to be incompressible, one must take 
care that the additional constraints do not become too restrict ive|§| Indeed, a simple 
argument shows that the model of an incompressible rod which rolls without sliding, 
and which cannot move transversally, can only move like a rigid body. 

There are two immediate solutions: either one relaxes the incompressibility constraint, 
or one allows the rod to move laterally as well. Either solution introduces a lot of 
mathematical tedium which greatly obscures the physical background of the system. 
For this paper, we will therefore consider a simplified model containing aspects of both 
models. 

In particular, we will assume that the motion of the nonholonomic rod is such that the 
incompressibility constraint is satisfied approximately throughout the motion; this is 
equivalent to the following assumption: 



vW + W=l. (23) 

This was pointed out to me by W. Tulczyjew and D. Zenkov. 
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By neglecting the incompressibility constraint in the Lagrangian, a simplified model 
is then obtained. Of course, this new model is a mathematical simplification of the 
true physics. However, numerical simulations show that \/{x') 2 + (y') 2 is bounded 
throughout the motion, and it's therefore reasonable that the dynamics of this model is 
close to the true dynamics. One could think of the mathematical model as describing a 
Cosserat rod whose constitutive equation is specified on mathematical grounds, rather 
than derived from first principles. 

The constraints of rolling without sliding are given by (see [251126] ): 

x + RO sin ip = and y - R6 cos <p = 0, (24) 

where R is the radius of the laminae. By eliminating the slope ip we then obtain 

x + My' = and y- R6x' = 0. (25) 

Incidentally, the passage from ( 1241) to (l25l) again illustrates why derivatives with respect 
to time play a fundamentally different role as opposed to the other derivatives. 

The Lagrangian density of the nonholonomic rod is still given by ( 1181) ; we recall that it 
is of second order, as the stored energy function (fl6l) is of grade two. The constraint 
on the other hand is of first order. By demanding that the action be stationary under 
variations compatible with the given constraint (a similar approach to section [273T) . we 
obtain the following field equations: 

Definition 3.2. A section <fi of it is a solution of the nonholonomic problem if and only 
iflmj 1 ^ C C, and, along C, 

(f<f>y(j 3 Vjn L ) = (26) 

for all tt -vertical vector fields VonY such that {j 1 <£>)*{ j l V_l a) = for all a G F. 



The left-hand side of (1261) is just the Euler-Lagrange equation ( TlOl) for a second-order 
Lagrangian. As the constraint is first order, it can be treated exactly as in section 12.31 
In coordinates, the nonholonomic field equations hence are given by 



dL d / dL \ d 2 f dL 



+ 



dy a dx^ ydy^J dxMx 1 ' ydy^ 



' hi/ 



(iV) = A a ^(jV) 



By substituting the Lagrangian (1T81) . without the inextensibility constraint, and the 
constraints ( 1231) into the Euler-Lagrange equations, we obtain the following set of 
nonholonomic field equations: 

px + Kx"" = X 

pij + Ky"" = n (27) 

a6-(36" = R{Xy' + fix') 

where A and /i are Lagrange multipliers associated with the nonholonomic constraints. 
These equations are to be supplemented by the constraint equations (1251) . 
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In the familiar case of the rolling disc, it is well known that energy is conserved. There 
is a similar conservation law for the nonholonomic rod. 

Proposition 3.3. The total energy ( fH|) is conserved for each solution of the 
nonholonomic field equations (F?\ ) and constraints W3\) . A fortiori, the solutions of 
the nonholonomic field equations satisfy the local conservation law d[(j 3 (p)*J[] = 0, 
where J[ is the momentum map associated to time translation introduced in (21\) . 



Proof: This follows immediately from proposition 16.11 in the appendix, and the fact 
that J^ (or rather its prolongation to J l ix) annihilates F along the constraint manifold. 
Indeed, the bundle of (n + l)-forms F is generated by Q 1 and $ 2 , defined as follows: 

$ x = (dx - xdt) A ds + Ry'(d6 - Odt) A ds; 
$ 2 = (dy - ydt) A ds - Rx'(d9 - 9dt) A ds. 

Therefore, we have that 



(§i) 1 J$1 = -( i + W)ds, 



which vanishes when restricted to C. A similar argument shows that the contraction of 
Jj: with $ 2 vanishes. Hence, proposition 16.11 can be applied; the associated momentum 
map is just (12T1) . □ 



4. Discrete nonholonomic field theories 

In this section we present an extension to the case of field theories of the discrete 
d'Alembert principle described in [27]. We also study an elementary numerical 
integration scheme aimed at integrating the field equations ( 1271) . 

As in the previous sections, we will consider the trivial bundle n with base space 
X = E x M (where M = [0, £]), and total space Y the product XxS, where S = R 2 x S 1 . 
Our discretization scheme is the most straightforward possible: The base space X will 
be discretized by means of the uniform mesh Z x Z, and the total space Y by replacing 
it with I x I 2 x 1. 

4-1. Discrete Lagrangian field theories 

We begin by giving an overview of discrete Lagrangian field theories, inspired by [211 128] . 
In order to discretize the second-order jet bundle, we need to approximate the derivatives 
of the field (of first and second order). This we do by means of central differences with 
spatial step k and time step h: 

V w ^ , V » 2k ' V W k 2 ' ^ ' 
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where r] stands for either x or y. Other derivatives will not be needed. For 8, we use 

« 9n+1 '\~ 6n ' 1 and e >n 9 ^+\- 9 ">\ (29) 



h k 

Let M. be the uniform mesh in X = R 2 whose elements are points with integer 
coordinates; i.e. M. = Z x Z. The elements of M. are denoted as (n, i), where the 
first component refers to time, and the second to the spatial coordinate. We define a 
9 -cell centered at (n, i) G /A, denoted by [x]( n ,i), to be a nine-tuple of the form 

[x\{n,i) '■= {{n - 1, * - 1), (n - 1, i), (n - 1, z); (n, i - 1), (30) 

(n, z), (n, z + 1); (n + 1, z - 1), (n + 1, i), (n + 1, i + 1)) 

It is clear from the finite difference approximations that a generic second-order jet j' 2 </> 
can be approximated by specifying the values of <fi at the nine points of a cell. 

However, in the case of the nonholonomic rod, the Lagrangian depends only on the 
derivative coordinates whose finite difference approximations were given in ( 1281) and 
(129]) . Therefore, we can simplify our exposition by defining a 6-cell at (n,i) to be the 
six-tuple 



[x]( n ,i) ■= ((n, z" - 1), (n, i), (n, z + 1); (n + 1, z - 1), (n + 1, i), (n + 1, z + 1)). 



(31) 



We will refer to 6-cells simply as cells. Let us denote the set of all cells by X 6 := {[x]( n) j) : 
(n, i) G M.}. We now define the discrete 2nd order jet bundle to be J|7r := X 6 x S 6 
(see [21| [28l [29]). A discrete section of 7r (also referred to as a discrete field) is a map 
: M. — > S*. Its second jet extension is the map j 2 : X f 



Jj7r defined as 



j 2 0(W(n,i)) : = (W(n,i); 0(»l), • • • 4>(Xq)), 



where xi,...,xq are the vertices that make up [#]( ni j) (ordered as in ( 1311) ). Given a 
vector field W on F, we define its second jet extension to be the vector field j 2 W on Jj 



given by 



fW([x]- A, . . . , / 6 ) = {W( Xl , h), W(x 2 , h 



Let us now assume that a discrete Lagrangian L d : Jjn - 
Sd is then defined as 

[x] 

Given a vertical vector field V on K and a discrete field 
family <p e by composing with the flow $ of V: 

£ ([x]) = ([x];$,(0([x]) 1 ),...,$ e (0([x]) 6 )). 



..,M/(x 6 ,/ 6 )). 

IR is given. The action sum 

(32) 
we obtain a one-parameter 

(33) 



The variational principle now consists of seeking discrete fields <fi that extremize the 
discrete action sum. The fact that 6 is an extremum of S under variations of the form 



A class of nonholonomic kinematic constraints in elasticity 22 

(I3"3"|) is expressed by 
J2 (X(<t>(n,0), DiL(j 2 (/»(N(n, i+ i))) + D 2 L(f<p([x} M )) + D 3 LO'V(N(n,<-i))) (34) 

(n,i)eM 

+ D 4 £0*V(M(n-i,i+i))) + D B L{j 2 <f>{[x] ir ^i,i))) + A$£(;V(M (n-i,i-i)))) = 0. 

As the variation X is completely arbitrary, we obtain the following set of discrete Euler- 
Lagrange field equations: 

D 1 L(f ( j ) ([x} inji+1) )) + D 2 L(f<j ) ([x} M )) + D 3 L(;'V(M(M-i))) + ( 35 ) 
D 4 L(j 2 (/.([x] (n _ lii+ i ) )) + I> 5 ^0*V(b](«-i,i))) + A>£(j 2 0(k](n-i,;-i))) = 0. 

for all (n, i). Here, we have denoted the values of the field (p at the points (n, i) as n> j. 

^.£. T/ie discrete d'Alembert principle 

Our discrete d'Alembert principle is nothing more than a suitable field-theoretic 
extension of the discrete Lagrange-d'Alembert principle described in [27]. Just as in 
that paper, in addition to the discrete Lagrangian L d , two additional ingredients are 
needed: a discrete constraint manifold Q C Jjff and a bundle of constraint forces Fd on 
Jj7r. However, as our constraints (in particular (125]) ) are not linear in the derivatives, 
as opposed to the case in [27], our analysis will be more involved. 

The discrete constraint manifold Cd "-^ J\k will usually be constructed from the 
continuous constraint manifold C by subjecting it to the same discretization as used 
for the discretization of the Lagrangian (i.e. (128]) and (129]) ). To construct the discrete 
counterpart Fd of the bundle of discrete constraint forces, somewhat more work is 
needed. 

Remark 4.1. For the discretization of the constraint manifold, it would appear that 
we need a discrete version of the first-order jet bundle as well. A similar procedure as 
for the discretization of the second-order jet bundle (using the same finite differences as 
in (128]) shows that a discrete 1-jet depends on the values of the field at the same four 
points of a cell as a discrete 2-jet: the difference between J\tx and Jjrr lies in the way 
in which the values of the field at these points are combined. Therefore, we can regard 
the discrete version of C, to be defined below, as a subset of Jjn. o 

4-2.1. The bundle of discrete constraint forces In this section, we will construct Fd 
by following a discrete version of the procedure used in section 12.31 Just as in the 
continuous case, it is here that the difference between spatial and time derivatives 
will become fundamental. Indeed, we will discretize with respect to space first, and 
(initially) not with respect to time. It should be noted that the construction outlined in 
this paragraph is not entirely rigorous but depends strongly on coordinate expressions. 
Presumably, one would need a sort of discrete Cauchy analysis in order to solidify these 
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arguments. For now, we will just accept that this reasoning provides us with the correct 
form of the constraint forces. 

For the sake of convenience, we suppose that C is given by the vanishing of k independent 
functions (p a on J 1 ^. By applying the spatial discretizations in ( l28i) and ( 1291) to <p a , 
we obtain k functions, denoted as ^1/2' on ^1 x TS. We define the semi-discretized 
constraint submanifold C1/2 to be the zero level set of the functions <fi/ 2 - 

Consider now the forms 

$? /2 := J*(d< /2 ) 

(where J is the vertical endomorphism on TS); they are the semi-discrete counterparts 
of the forms $ a defined in (1151) . The forms $?, 2 are semi-basic. By discretizing the time 
derivatives, however, we obtain a set of basic forms on J%tt, which we also denote by 
$"/ 2 . An example will make this clearer. 

Example 4.2. Consider, for instance, the constraint manifold C ^-> J 2 n defined as the 
zero level set of the function if = A^y^y^ + B b (y\) 2 , where A ab and B b are constants. 
By applying (1281) and (|29|) , it follows that Cd is given as the zero level set of the function 

Miv}) ■= A ab - - + B b y — 

for [y] G J^7r, and C\/ 2 as the zero level set of the function 

/r 1 \ a -a Vn,i+1 Vn,i-1 ,-> / Vn,i+1 Vn,i-1 

<Pi/2{\y],v) :=A ab v — + B b \ — 

for [y] G J^7r and v G TS. The bundle Fj is then generated by the one-form 
$1/2 := S*(d(p 1/2 ), or explicitly, 

$ = 4 /"'* 1 ~ y " ,i - 1 dy a . 

4-2.2. The discrete nonholonomic field equations Assuming that Ld, Cd and Fd are 
given (their construction will be treated in more detail in the next section), the derivation 
of the discrete nonholonomic field equations is similar to the continuum derivation: we 
are looking for a discrete field <fi such that Imj 1 ^ C Cd and such that (f) is an extremum 
of ( 132]) for all variations compatible with the constraints, in the sense that the variation 
X satisfies, for all (n,i), 

From fl34l) we then obtain the discrete nonholonomic field equations: 
L>iL(j 2 0(M(n,*+i))) + D 2 L(f(j)([x] (n , i) )) + D 3 LO' 2 0(M(n,i-i))) + D A L(J 2 (f>([x]( n -i t i+i))) 
+ D 6 £0'V(M(n-i,i))) + As£0'V(M(»-i,i-i))) = A a $? /2 (j 2 0(M(n,i))), (36) 
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where the Lagrange multipliers X a are to be determined from the requirement that 
lmj 2 C C d . 



4-3. An explicit, second-order algorithm, 

In this section, we briefly present some numerical insights into the nonholonomic field 
equations of section 13.51 Our aim is twofold: for generic boundary conditions, the 
nonholonomic field equations ( 1271) cannot be solved analytically and in order to gain 
insight into the behaviour of our model, we therefore turn to numerical methods. 
Secondly, in line with the fundamental tenets of geometric integration, we wish to show 
that the construction of practical integration schemes is strongly guided by geometric 
principles. 

In discretizing our rod model, we effectively replace the continuous rod by N rigid rolling 
discs interconnected by some potential (see [30]). This is again an illustration of the fact 
that the constraints are truly nonholonomic. Our integrator is just a concatenation of 
the leapfrog algorithm for the spatial part, and a nonholonomic mechanical integrator 
for the integration in time. 

As a first attempt at integrating (1271) . we present an explicit, second-order algorithm. 
In the Lagrangian, the derivatives are approximated by 

x ~ : and x ~ — , 

h k 2 

where h is the time step, and k is the space step. Similar approximations are used for 
the derivatives of y, and for 9 we use 

h ~, "n+l,i "n,i j n t _ "n,i+l "n,i /o>7\ 

w — and « — ' — -. (37) 

h k 

The discrete Lagrangian density can then be found by substituting these approximations 
into the continuum Lagrangian (fi~8l) . Explicitly, it is given by 

d = ~Oh2 \\ Xn + 1 < i _ X n,i) + \Vn+l,i ~ Dn,i) ) + ^T^{"n+l,i ~ Vn,i) ~ ^T^{^n,i+1 ~ Un,i) 



K K 

— 77T^(^n,i-l — 2x n) j + X n j+i) — — z-^\y n ,i-l ~ ^Vn,i + 2/n,i+l) • (38) 

Note that La only depends on four of the six points in each cell (see (}3~Tj) ). The discrete 
constraint manifold Cd is found by discretizing the constraint equations (125]) . In order 
to obtain a second-order accurate approximation, we use central differences: 

/ -EnA+l %n,i— 1 

rp i~*~l : ; 

2k 
(and similar for y', x, y, 6) and hence we obtain that Cd is given by 

%n+l,i — ^n-l,i + ^r(^n+l,i _ &n-l,i) (Z/n,i+l — Un,i-l) = 0, (39) 
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and 

Vn+l,i ~ Vn-l,i ~ 7Tr(#n+M ~ $n-l,i) \%n t i+l — #n,i-l) = 0, (40) 

for all (n,i). The semi-discrete constraint manifold C1/2, on the other hand, is given by 

p 

•Kn,i * n j, w,i \Vn,i-\-l i/n,i— 1/ U, 



and 



p 

Vn/i TTT Vn,i\%n,i+1 Xn,i—X) U, 



and hence Fd is generated by 



p p 

$! = dx + — (y n ,j + i - Z/rM-i)d0 and $ 2 = dy - — (i n , i+ i - ar n< _i)d0. 
2A; 2k 



We conclude that the discrete nonholonomic field equations (|36l) are in this case 

h 2 K 
x n+1>i - 2x n ,i + x„_ M = /i 2 A, - yA 4 i n)j (41) 



and 

as well as 



y n+1;i - 2y n>i + y n _i 4 = /i 2 /ii n-A 4 y n>i (42) 



//) o/Q i a \ d(,2 / \ Un,i+1 Un,i— 1 x n,i+l x n,i—l \ . P* 1 A 2/i 

a^n+i^ — /f7 nj j + V n -ij) = nil Ai — /ij — + -j-^-A f ni j, 



2A; ^ 2fc / fc 2 

where A 2 and A 4 are the 2nd and 4th order finite difference operators in the spatial 
direction, respectively: 

^ Jn.i ■ Jn,i+1 ^Jn,i ~r Jn,i—1 

and 

A f n .i := J n ,i+2 ~ 4/ n ,i+l + 6/ nj j — 4/ n j__ 1 + Jn : i-2- 

In order to determine Aj and //», these equations need to be supplemented by the discrete 
constraints (1391) and (1401) . 

For the purpose of numerical simulation, the following values were used: a = 1, (3 = 0.8, 
p = 1, K = 0.7, £ = 4, and i? = 1. For the spatial discretization, 32 points were used 
(corresponding to k « 0.1290) and the time step was set to h = l/8k 2 , a fraction of the 
maximal allowable time step for the Euler-Bernoulli beam equation (see [31 j ) . The ends 
of the rod were left free and the following initial conditions were used: 

71 7TS 

r o (s) = {s,0), 9 (s) = --cos — and r (s) = (0,0), 6 (s) = 0. 
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Figure 3. Motion of the rod from t = to t « 4.5. 



An mpeg movie (created with Povray, an open source raytracer) depicting the motion 
of the nonholonomic Cosserat rod is available from the author's web pagqjjj. In figure [3j 
an impression of the motion of the rod is given. The arrows represent the director field 
d 3 and serve as an indication of the torsion. The rod starts from an initially straight, 
but twisted state and gradually untwists, meanwhile effecting a rotation. 

In figure HJ the energy of the nonholonomic rod is plotted. Even though our algorithm 
is by its very nature not symplectic (or multi-symplectic - see [32]), there is still the 
similar behaviour of "almost" energy conservation. 



http: //users .ugent . be/$\sim$jvkersch/nonholonomic/ 
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Figure 4. Energy behaviour of the integration algorithm on moderate time scales (the 
interval [0,150]). 



5. Conclusions 

It is clear that the study of nonholonomic field theories forms a vast subject. This 
paper gives only a brief survey of a number of straightforward results, but there are 
many more things to be explored. An acute problem is the lack of an extensive number 
of interesting examples; while this of course impedes progress on the theoretical front, 
there are nevertheless a number of points worth investigating, which we now discuss. 

In proposition 13.31 we used the fact that the bundle of reaction forces is annihilated by 
the generator of time translations in order to prove conservation of energy. Even when 
this is not the case, experience from mechanics (see [33" l l3"4" l l35]) as well as from different 
types of nonholonomic field theories (see [36]) seems to suggest that it might be possible 
to prove a nonholonomic momentum equation instead. 

From a numerical point of view, the explicit algorithm of section 14.31 is not very 
accurate. It is second-order in space and time and suffers from a restrictive stability 
condition. The development of more sophisticated integration schemes that exactly 
preserve the nonholonomic constraints would definitely be very interesting. Perhaps 
the most interesting of all, at least in line with the current investigations, would be a 
simulation of a nonholonomic model with a more physical constitutive equation than 
the one used in section 13.51 

6. Appendix: the nonholonomic Noether theorem 



In the derivation of the Euler-Lagrange equations, both in the free as in the constrained 
case, we have restricted our attention to vertical variations. While it is well known 
that arbitrary variations do not yield any new information beyond the Euler-Lagrange 
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equations (see [5]), the situation is not at all clear for nonholonomic field theories. 

This is especially important for the derivation of the nonholonomic Noether theorem. 
Therefore, we propose the following modified bundle of reaction forces: if F — (A^dy "), 
then 

F := {A a a (dy a - y^) A d n x ) C ^ +1 )(jV). (43) 

This situation is reminiscent of the comparison between Bridges' (n+1) multisymplectic 
1-forms and the Cartan (n + l)-form 0^ in the work of Marsden and Shkoller [37], and 
is similar to the variational derivation of the Cartan form: if only vertical variations are 
taken into account, then the Cartan form is missing the d n+1 x term (see [3]). 

The precise form of the bundle F can be derived by using arguments from Cauchy 
analysis (see [IE]). The idea is to reformulate the field equations as a mechanical 
system on an infinite dimensional configuration space. The reaction forces can then be 
introduced in a straightforward way on this infinite-dimensional space, and by returning 
to the jet bundle one then obtains ( 143}) . The details of that derivation would lead us too 
far; more information on this technique will appear in a forthcoming publication (see 
also |38j). For now, we will simply accept the bundle F as given. 

In [36] a similar type of bundle was used in the derivation of the nonholonomic 
momentum lemma. From that paper, we cite the following nonholonomic Noether 
theorem: 

Proposition 6.1. Let L be a G-invariant Lagrangian density. Assume that £ G g is 
such that £ji,rJ a = along C for all a G F . Then the following conservation law holds: 

d[0'V) *4\ = o, 

for all sections <f> of n that are solutions of the nonholonomic field equations. 

Note that a vertical vector v belongs to F° if and only if (v, a) = for all a G F. 
Only for non-vertical vectors there is a difference between F and F. Therefore, if £ji n 
in proposition 16. II is vertical, then the nonholonomic Noether theorem follows from the 
Euler-Lagrange equations (j2.9p . In the other case, the techniques from [36] have to be 
used. 
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